﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using NGDTConsole.Data.Enumerations;

namespace NGDTConsole.Data.Parameters
{
    public class Ellipsoid
    {        
        private double refLatitude = 0.0;
        private double refLongitude = 0.0;
        private double parametricLatitude = 0.0; //u
        private double majorAxis = 0.0;
        private double minorAxis = 0.0;
        
        private double ellipsoidEccentricity; //e
        private double refSecEccentricity; //ePrime
        private double ellipsoidFlattening; //f        
        private double radiusOfCurvaturePrimeVertical; //P
        private double radiusOfCurvatureMeridian;


        public Ellipsoid(GeodeticMeasurementUnits unit, double majorAxis, double minorAxis, int precision = 0)
        {
            //assert the values
            bool n1 = Double.IsNaN(refLatitude) || Double.IsInfinity(refLatitude);
            bool n2 = Double.IsNaN(refLongitude) || Double.IsInfinity(refLongitude);
            bool n3 = Double.IsNaN(majorAxis) || Double.IsInfinity(majorAxis);
            bool n4 = Double.IsNaN(minorAxis) || Double.IsInfinity(minorAxis);

            
            if (!n1 && !n2 && !n3 && !n4)
            {                
                this.majorAxis = Math.Round(majorAxis, precision);

                //set the primary variable for other parameters                
                this.ellipsoidFlattening = CalculateFlattening(precision);
                
                //set other derived variables                
                this.ellipsoidEccentricity = CalculateEccentricity(precision);

                this.refSecEccentricity = CalculateEllipsoidSecEccentricity(precision);
            }
        }

        public Ellipsoid(GeodeticMeasurementUnits unit, double majorAxis, double flattening, double refLatitude = 0, double refLongitude = 0, int precision = 0)
        {
            bool n1 = Double.IsNaN(refLatitude) || Double.IsInfinity(refLatitude);
            bool n2 = Double.IsNaN(refLongitude) || Double.IsInfinity(refLongitude);
            bool n3 = Double.IsNaN(majorAxis) || Double.IsInfinity(majorAxis);

            if (unit.Equals(GeodeticMeasurementUnits.Radian))
            {
                this.refLatitude = refLatitude;
                this.refLongitude = refLongitude;
            }
            if (unit.Equals(GeodeticMeasurementUnits.DecimalDegree))
            {
                this.refLatitude = Math.Round((refLatitude * Math.PI / 180), precision);
                this.refLongitude = Math.Round((refLongitude * Math.PI / 180), precision);
            }
            if (!n1 && !n2 && !n3)
            {                
                this.majorAxis = Math.Round(majorAxis, precision);
                this.ellipsoidFlattening = Math.Round(flattening, precision);
                
                //calculate the derived values in right order
                this.ellipsoidEccentricity = CalculateEccentricity(precision);
                this.refSecEccentricity = CalculateEllipsoidSecEccentricity(precision);
                this.radiusOfCurvaturePrimeVertical = CalculatePrimeVertical(this.refLatitude, precision);
                this.radiusOfCurvatureMeridian = CalculateCurvatureMeridian(this.refLatitude, precision);
            }
        }

        #region Ellipsoid Properties
        
        public double MajorAxis
        {
            get { return this.majorAxis; }
            set { this.majorAxis = value; }
        }

        public double MinorAxis
        {
            get { return this.minorAxis; }
            set { this.minorAxis = value; }
        }

        public double Eccentricity
        {
            get
            {
                return CalculateEccentricity();
            }
        }

        public double SecondEccentricity
        {
            get
            {
                return CalculateEllipsoidSecEccentricity();
            }
        }

        public double ParametricLatitude
        {
            get { return this.parametricLatitude; }
            set { this.parametricLatitude = value; }
        }
        /*
         * a is the ellipsoidal semi-major axis
b is the ellipsoidal semi-minor axis
f is the flattening of the ellipsoid where 1/f = a/(a - b)
        
         */
        //derived properties
        public double Flattening
        {
            get
            {
                return this.ellipsoidFlattening;
            }
            set
            {
                this.ellipsoidFlattening = value;
            }
        }

        public double InverseFlattening
        {
            get
            {
                return 1 / this.ellipsoidFlattening;
            }
        }
       
        /*
         * ρ is the radius of curvature of the meridian at latitude ϕ,
where ρ = a(1 – e2)/(1 – e2sin2ϕ)3/2
ν is the radius of curvature on the prime vertical (i.e. perpendicular to the meridian) at
latitude ϕ, 
         * 
         */
        public double MeridianCurvature
        {
            get
            {
                return this.radiusOfCurvatureMeridian;
            }            
        }

        public double PrimeVertical
        {
            get
            {
                return this.radiusOfCurvaturePrimeVertical;
            }
            set { this.radiusOfCurvaturePrimeVertical = value; }
        }

        #endregion

        #region Private method of Ellipsoid

        private double CalculateFlattening(int fraction = 0)
        {
            //assert the values
            bool n1 = Double.IsInfinity(this.minorAxis) || Double.IsInfinity(this.majorAxis);
            bool n2 = Double.IsNaN(this.majorAxis) || Double.IsNaN(this.minorAxis);
            if (!n1 && !n2)
            {
                double cal = (this.majorAxis - this.minorAxis) / this.majorAxis;
                if (fraction > 0)
                {
                    cal = Math.Round(cal, fraction);
                }
                return cal;
            }
            else
            {
                return 0.0;
            }
        }

        //The eccentricity is defined as eccentricitySqauare
        private double CalculateEccentricity(int fraction = 0)
        {
            if (this.ellipsoidFlattening > 0 || !Double.IsNaN(this.ellipsoidFlattening))
            {
                double cal = this.ellipsoidFlattening * (2 - this.ellipsoidFlattening);
                if (fraction > 0)
                {
                    cal = Math.Round(cal, fraction);
                }
                return cal;
                
            }
            else
            {
                return 0.0;
            }
        }

        //e' is the second eccentricity where e'2 = e2 /(1 –e2)
        private double CalculateEllipsoidSecEccentricity(int fraction = 0)
        {
            bool n1 = (!Double.IsNaN(this.ellipsoidFlattening)) && (!Double.IsInfinity(this.ellipsoidFlattening));
            if (n1)
            {                
                double ePrimeSqr = (this.ellipsoidFlattening * (2 - this.ellipsoidFlattening)) / Math.Pow((1 - this.ellipsoidFlattening), 2);
                if (fraction > 0)
                {
                    ePrimeSqr = Math.Round(ePrimeSqr, fraction);
                }
                return ePrimeSqr;
            }
            else
            {
                return 0.0;
            }
            
        }
        
        private double CalculateCurvatureMeridian(double refLatitude, int fraction = 0)
        {
            if (Eccentricity > 0 || !Double.IsNaN(Eccentricity))
            {
                
                double e = Eccentricity;
                double a = MajorAxis;
                double num = a * Math.Pow((1 - e), 2);
                double denum = Math.Pow(1 - (e * (Math.Pow(Math.Sin(refLatitude), 2))), 1.5);
                double m = num / denum;
                if (fraction > 0)
                {
                    m = Math.Round(m, fraction);
                }
                return m;
            }
            else
            {
                return 0.0;
            }
        }
        //where ν = a /(1 – e2sin2ϕ)1/2
        private double CalculatePrimeVertical(double refLatitude, int fraction = 0)
        {
            if (!Double.IsInfinity(Eccentricity) && !Double.IsNaN(Eccentricity))
            {
                double a = MajorAxis;
                double e = Eccentricity;
                double v = 0.0;
                v = a / (Math.Sqrt(1 - (e * (Math.Pow(Math.Sin(refLatitude), 2)))));
                if (fraction > 0)
                {
                    v = Math.Round(v, fraction);
                }
                return v;
            }
            else
            {
                return 0.0;
            }
        }

        #endregion
    }
}
